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Abstract We propose a new image denoising algorithm when the data is con- 
taminated by a Poisson noise. As in the Non-Local Means filter, the proposed 
algorithm is based on a weighted linear combination of the observed image. 
But in contract to the latter where the weights are defined by a Gaussian ker- 
nel, we propose to choose them in an optimal way. First some " oracle" weights 

■ are defined by minimizing a very tight upper bound of the Mean Square Error. 
For a practical application the weights are estimated from the observed image. 
We prove that the proposed filter converges at the usual optimal rate to the 

, true image. Simulation results are presented to compare the performance of 

■ the presented filter with conventional filtering methods. 
A-^ ■ 

I Keywords Poisson noise • Mean Square Error • oracle estimate • Optimal 

C/3 ■ Weights Filter 



^ ■ 1 Introduction 

In a variety of applications, ranging from nuclear medicine to night vision and 
, from astronomy to traffic analysis, data are collected by counting a series of 

discrete events, such as photons hitting a detector or vehicles passing a sensor. 
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Many such problems can be viewed as the recovery of the intensity from the 
indirect Poisson data. The measurements are often inherently noisy due to 
low count levels, and we wish to reconstruct salient features of the underlying 
phenomenon from these noisy measurements as accurately as possible. 

There are many types of methods to reconstruct the image contaminated by 
the Poisson noise. The most popular method is performed through a Variance 
Stabilizing Transform (VST) with the following three-step procedure. First, 
the variance of the Poisson distribution is stabilized by applying a VST. So that 
the transformed data are approximately homoscedastic and Gaussian. The 
VST can be an Anscombe root transformation (Anscombe Q and Borovkov 
0), multiscal VSTs (Bardsley and Luttman [33,]), Conditional Variance Sta- 
bilization (CVS) (Jansen [31), or Haar-Fisz transformation (Fryzlewicz and 
Nason [13, ini). Second, the noise is removed using a conventional denois- 
ing algorithm for additive Gaussian white noise, see for example Buades, Coll 
and Morel (2005 [3), Kervrann (2006 [13]), Aharon and Elad and Bruck- 
stein (2006 [i|), Hammond and Simoncelh (2008 [13]), Polzehl and Spokoiny 
(2006 [^), Hirakawa and Parks (2006 [i^l), Mairal, Samro and Elad (2008 
[13), Portilla, Strela, Wainwright and Simoncelh (2003 [23)^oth and Black 
(2009 [ii]), Katkovnik, Foi, Egiazarian, and Astola (2010 Dabov, Foi, 

Katkovnik and Egiazarian (2006 [^), Abraham, Abraham, Desolneux and Li- 
Thiao-Te (2007 [1,]), and Jin, Grama and Liu (2011 [H). Third, an inverse 
transformation is applied to the denoised signal, obtaining the estimate of the 
signal of interest. Makitalo and Foi (2009 |2lj and 2011 122[) focus on this last 
step, and introduce the Exact Unbiased Inverse (EUI) approach. Zhang, Fadili, 
and Starck (2008 [33|]), Lefkimmiatis, Maragos, and Papandreou (2009 (isj), 
Luisier, Vonesch, Blu and Unser (2010 [l^l) improved both the stabilization 
and the inverse transformation. 

Regularization based on a total variation scminorm has also attracted sig- 
nificant attention, see for example Beck and TebouUe (2009 [^), Bardsley and 
Luttman (2009ji|), Setzer, Steidl and Teuber (2010 [3l|). Nowak and Ko- 
laczyk (1998 [24| and 2000 [2^) have investigated reconstruction algorithms 
specifically designed for the Poisson noise without the need of VSTs. 

In this paper, we introduce a new algorithm to restore the Poisson noise 
without using VST's. We combine the special properties of the Poisson dis- 
tribution and the idea of Optimal Weights Filter [l5|| for removing efficiently 
the Poisson noise. The use of the proposed filter is justified both from the 
theoretical point of view by convergence theorems, and by simulations which 
show that the filter is very effective. 

The paper is organized as follows. Our main results are presented in Section 
[5] where we construct an adaptive estimator and give an estimation of its rate 
of convergence. In Section |31 we present our simulation results with a brief 
analysis. Proofs of the main results are deferred to Section HI 
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2 Construction of the estimator and its convergence 

2.1 The model and the notations 

We suppose that the original image of the object being photographed is a 
integrable two-dimensional function f{x), x G (0,1] x (0,1]. Let the mean 
value of / in a set B^; be 

A{B,) ^N^ J J{t)dt. 

Typically we observe a discrete data set of counts Y — {7V(Ba;)}, where A/'(Bx) 
is a Poisson random variable of intensity yl(Bx). We consider that if B^^nBj^ = 
0, then NiBx) is independent of NiBy). For a positive integer TV the uniform 
N X N grid on the unit square is defined by 

r 1 2 TV- 1 1 ^ 
^^b'TV'--'^'^} ■ 

Each element x of the grid I is called pixel. The number of pixels is n = N'^. 
Suppose that x = {x^^\x^^^) £ I, and B^^ = (a;*^) - 1/N,x^^')] x (a;*^) - 
1/A^, x'^^^]. Then {Ba;}^:^! is a partition of the square (0, 1] x (0, 1]. The image 
function / is considered to be constant on each B^;, x £ 1. Hence we get a 
discrete function f{x) — ACBx), x € I. The denoising aims at estimating the 
underlying intensity profile f{x). In the sequence we shall use the following 
important property of the Poisson distribution: 

E(AA(B,)) - Var(AA(B,)) = f{x). (2) 

Actually the Poisson noise model can be viewed as the following additive noise 
model 

Y{x)^f{x)+e{x), (3) 

where 

eiy) = Y{x) - fix). (4) 

may be considered as an additive heteroscedastic noise related to the Poisson 
model. Due to ©, we have E(e(y)) = and Yar{e{y)) Yar{Y{y)) = f{x). 
Let us set some notations to be used throughout the paper. The Euclidean 

norm of a vector x = (xi, Xd) G R"* is denoted by ||x||2 = (SiLi ) ^ • The 
supremum norm of x is denoted by | 

•^lloo — ^^Pi<i<d l*^*! ■ The cardinality of 
a set A is denoted card A. For any pixel xq £ \ and a given /i > 0, the square 
window 

'^xo,h^ {x £l:\\x - xoWoo <h} (5) 

is called search window at xq. We naturally take ft. as a multiple of (h = 
for some fc e {1, 2, • • • , N}). The size of the square search window \Jxo,h is the 
positive integer number 

M = {2Nh + 1)2 = card U^„j,. 
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For any pixel x e ^xo,h and a given 77 > 0. Consider a second square window 
\Jx.ri of size 

m = {2Nr] + 1)^ = card Uxo,n- 

We shall call U^;,,, local patches and XJx.h search windows. Finally, the positive 
part of a real number a is denoted by : 

^ ( a a a > 0, 
" \ if a < 0. 



2.2 Construction of the estimator 

Let /i > be fixed. For any pixel G I consider a family of weighted estimates 
fh,wixo) of the form 

fh.w{xo) = ^ w{x)Y{x), (6) 

where the unknown weights satisfy 

w(x) > and w{x) = 1. (7) 

The usual bias and variance decomposition of the Mean Square Error gives 
E (^h.^xo) - f{xo)y - Bias' + Var, (8) 

where 
and 

Var — E^ w(a;)^/(a;). 

The decomposition ([8]) is commonly used to construct asymptotically minimax 
estimators over some given classes of functions in the nonparametric function 
estimation. With our approach the bias term Bias^ will be bounded in terms 
of the unknown function / itself. As a result we obtain some " oracle" weights 
w adapted to the unknown function / at hand, which will be estimated further 
using data patches from the image Y. 

First, we shall address the problem of determining the "oracle" weights. 
With this aim denote 



Pix) = Pf.xo (x) = \f{x) ~ f{xo)\ 



(9) 
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Note that the value of Pf.xo {x) characterizes the variation of the image bright- 
ness of the pixel x with respect to the pixel xq. From the decomposition ([5]), 
we easily obtain a tight upper bound in terms of the vector Pf^xo ■ 

E (A(xo) - fixo))' < gp,,^,H = gp{w), (10) 

where 

9pM=( E ^i^)p(^n + E ^i^)^fi^)- (11) 

From the following theorem we can obtain the form of the weights w which 
minimize the function gp{w) under the constraints ([7]) in terms of p {x) . For the 
sake of generality, we shall formulate the result for an arbitrary non-negative 
function p(x), x G Ux.h, not necessarily defined by 

Introduce into consideration the strictly increasing function 

^''W= E -7^Mx){t-p{x))+, t>0. (12) 

Let Kti- be the usual triangular kernel: 

i^tr (i) = (1 - |t|)+ , teR\ (13) 

Theorem 1 Let p{x) , x G \^xo,h be an arbitrary similarity function and let 
gp(w) be given by ill]) . Suppose that f{x) > for all x £ Vxo.h- Then there 
are unique weights which minimize gp{w) subject to Q), given by 

Ktr(^)/f{x) 

E,eu.„,.^*.(^)//(y) 

where a > Q is the unique solution of the equation 

Mp {a) - 1. (15) 
The proof of Theorem [T] is deferred to Section 14.11 
Remark 1 The bandwidth a > is the solution of 

E 7^P(^)(a- P(2^))^ = 1' 

and can be calculated as follows. We sort the set {p{x) \ x G ^^xaJi} in the 
ascending order — pi < p2 < • • • < Pm < Pm+i = +oo, where M — 
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card \Jxo,h- Let fi be the corresponding value of f{x) (we have fi = f{x) if 
Pi = p{x), X e l^xo.h)- Let 

ak = ^ , l<k<M, (16) 

and 

k* — maxjl < k < M \ ak > Pk] 

= niin{l < fc < M I flfe < p^} - 1, (17) 

with the convention that ak = oo ii pk — Q and that niin = M + 1 . The 
bandwidth a > can be expressed as a = a^. . Moreover, k* is also the unique 
integer k e {1, • • • , Af} such that a^ > Pfe and afe+i < Pk+i if fc < M . 

The proof of Remark [1] can be found in [l^. 

Let p(x) , X £ '^xa.h, be an arbitrary non-negative function and let Wp be 
the optimal weights given by (ITU) . Using these weights Wp we define the family 
of estimates 

//:(^o)= E ^p(^)^(^) (18) 

depending on the unknown function p. The next theorem shows that one can 
pick up an useful estimate from the family if the function p is close to the 
"true" function Pf,xo{x) = \f {x) - f {xa)\ , i.e. if 

p(x) = |/(x)-/(xo)|+<5„, (19) 

where i5„ > is a small deterministic error. We shall prove the convergence of 
the estimate under the local Holder condition 

\f{^) - f{y)\ <L\\x~ y\\L Vx, y G Vx^h+r,, (20) 

where /3 > is a constant, h > 0, -q > and xq S I. 

In the following, Ci > (z > 1) denotes a positive constant, and 0(a„) 
{n > 1) denotes a sequence bounded by c • a„ for some constant c > and 
all n > 1. All the constants Cj > and c > depend only on L and /?; their 
values can be different from line to line. Let 

r > ma.x{f{x) -.xel} (21) 

be an upper bound of the image /. 

Theorem 2 Assume that h > co?i~" with < a < 2/3+2 '^^^ ''O ^' 

that h = Con" 2/3+2 with cq > ci = (^F ^^"'"§^^2^"'"^'* ^ ''^ ■ Suppose also that 
the function / > satisfies the local Holder condition Let f^ixo) be 
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given by U8\). where the weights Wp are defined by ^4\ l <^^d, U5\} with p [x) — 
1/ {x) - f {xq)\ + Sn and dn^O (n^'^'j . Then 

EiKM-fMf = 0{n-^y (22) 



For the proof of this theorem see Section [ 

_ 1 

Recah that the bandwidth h of order n 2+2? is required to have the opti- 
mal minimax rate of convergence O (^n~ j of the Mean Square Error for 

estimating the function / of local Holder smoothness /3 (cf. e.g. 9]). To bet- 
ter understand the adaptivity property of the oracle fiixo), assume that the 
image / at xq has local Holder smoothness (3 (see [32|) and that h > cqtt."" 

with < a < 2/3+2 ' '^liich means that the radius ft, > of the search window 

_ 1 

UxQ,h has been chosen larger than the "standard" n ^3+2 , Then, by Theorem 
[51 the rate of convergence of the oracle is still of order n^^+^ . If we choose a 
sufhciently large search window Uxgjn then the oracle f^i^o) will have a rate 
of convergence which depends only on the unknown maximal local smoothness 
/? of the image /. In particular, if /3 is very large, then the rate will be close to 
rt~^/^, which ensures a good estimation of the flat regions in cases where the 
regions are indeed flat. More generally, since Theorem [5] is valid for arbitrary 
/?, it applies for the maximal local Holder smoothness at xq, therefore the 



oracle f^{xo) will exhibit the best rate of convergence of order n ^+^'^^0 at xq. 
In other words, the procedure adapts to the best rate of convergence at each 
point xo of the image. 

We justify by simulation results that the difference between the oracle 
computed with p{x) = pf^xoi^) = |/(a;) — f {^q)\ i and the true image /, is 
extremely small (see Table [T]). This shows that, at least from the practical 
point of view, it is justified to optimize the upper bound gpj, (w) instead of 
optimizing the Mean Square Error E {/^{xq) — f{xo))^ itself. 

The estimate with the choice p {x) — Pf\xo (x) will be called oracle filter. 
In particular for the oracle filter , under the conditions of Theorem [51 we 
have 

Eir^ixo)- f{xo)f <gp {wp)<cn-^. 

Now, we turn to the study of the convergence of the Optimal Weights Filter. 
Due to the difficulty in dealing with the dependence of the weights we shall 
consider a slightly modified version of the proposed algorithm: we divide the 
set of pixels into two disjoint parts, so that the weights are constructed from 
one part, and the estimation of the target function is a weighted mean along 
the other part. More precisely, we proceed as follows. Assume that xq S I. 
Denote 

1^0 - ja^o + (^-^, e I : I + i is pair 

and = Denote U^^^^ = Vx„h n I^^, and U^^,, = U,,, n V^^. Since 

E|y(a;) — Y{xo)\'^ = \f{x) — /(a;o)P + fixo) + fix), an obvious estimate of 
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E \Y{x) — Y{xo)\'^ is given by 



\— J2 \yiy)-Y{Ty)\\ 



cardU-.„„ 



where T — T^^.x is the translation mapping: Ty = x + {y — xq). Define an 
estimated similarity function p^g by 



^s^dW^ S iyM-y(Tv)A -v^) . (23) 



Pxo (x) = 



where 



The Optimal Weights Poisson Noise Filter (OWPNF) proposed in this paper 
is defined by 

hixo)^ w{x)Y{x), (24) 



^0. 



where 



argmin ^ w{x)pxo{x, xq) \ + fixo) ^ u;^(x). (25) 



In the next theorem, we prove that with the choice h — cqu 23+2 and 
r] = C2n^ , the Mean Square Error of the estimator fh{xo) converges nearly 

2/3 

at the rate n which is the usual optimal rate of convergence for a given 

Holder smoothness /3 > (see e.g. Fan and Gijbels (1996 [§|)). 

1 

Theorem 3 Assume that h = cqu'^^ with cq > ci = (ri^±|^||±2i^ ^ 

and that 7] = C2n^^3+2. Suppose that the function f{x) > satisfies the local 
Holder condition 12(A). Then 



nlhixo) - f{xo)f = O (n-^ In^ n) . (26) 
For the proof of this theorem see Section 14.31 
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3 Simulation 

For simulations we use the following usual set of 256x256 images: Spots[0.08, 4.99], 

Galaxy[0, 5], Ridges[0.05, 0.85], Barbara[0.93, 15.73] and Cells [0.53, 16.93] (see 

the first row of Figure[T|). All the images are included in the package "Dcnoising 

software for Poisson data" which can be downloaded at http://www.es. tut. fi/^foi/invansc/ 

We first do the simulations with the oracle filter which shows excellent visual 

quality of the reconstructed image. We next present our denoising algorithm 

and the numerical results which are compared with related recent works ( (22j 



and [33|). Each of the aforementioned articles proposes an algorithm specif- 
ically designed for Poisson noise removal (EUI+BM3D, MS- VST + 7/9 and 
MS- VST + B3 respectively). 

We evaluate the performance of a denoising filter / by using the Normalized 
Mean Integrated Square Error (NMISE) defined by 



NMISE : 




(/(^) - fix))' 
fix) 



where f{x) are the estimated intensities, f{x) are the respective true vales, 
and n* = card{/(a;) : f{x) > 0, a; £ I}. 



3.1 Oracle Filter 

In this section we present the denoising algorithm called Oracle Filter, and 
show its performance on some test images. 

Algorithm: Oracle Filter 

Repeat for each G I 

Let a = 1 (give the initial value of a) 
compute pixi) by ([27]) 
reorder pixi) as increasing sequence 
loop from fc = 1 to M 

else quit loop 
else continue loop 
end loop 

{a-p(x,)) + /f{x,) 



compute wix,) = (a-p(.,))+//(.,) 
compute /^(xo) = Ea;,GU,„ ^ w{xi)Y{x^). 



We calculate the optimal weights from the original image and compute the 
oracle estimate from the observed image contaminated by the Poisson noise. 
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Table 1 NMISE values when oracle estimator /j* is applied with different values of M. 



yize 


7x7 


9x9 


11 X 11 


13 X 13 


15 X 15 


17 X 17 


19 X 19 


Spots[0.08,4.99J 


0.0302 


0.0197 


0.0166 


0.0139 


0.0112 


0.0098 


0.0104 


Galaxy[0,5] 


0.0284 


0.0208 


0.0165 


0.0144 


0.0122 


0.0107 


0.0093 


Ridges[0.05,0.85] 


0.0239 


0.0178 


0.0131 


0.0109 


0.0098 


0.0085 


0.0074 


Barbara[0.93,15.73] 


0.0510 


0.0399 


0.0304 


0.0248 


0.0208 


0.0195 


0.0174 


CellsfO.53,16.931 


0.0422 


0.0323 


0.0257 


0.0216 


0.0191 


0.0164 


0.0146 




(a) Spots (b) Galaxy (c) Ridges (d) Barbara (e) Cells 



Fig. 1 The first row is the test images original, and the second row is the images restored 
by Oracle Filter with M = 19 X 19. 



For choosing the convenient size of the search windows, we do numerical ex- 
periments with different window sizes (see Table [ij. The results show that the 
difference between the oracle estimator and the true value / is extremely 
small. In Figure [TJ the second row illustrates the visual quality of the restored 
images by the Oracle Filter with M = 19 x 19. We can see that almost all the 
details have been retained. 



3.2 Performance of the Optimal Weights Poisson Noise Filter 

Throughout the simulations, we use the following algorithm for computing 
the Optimal Weights Poisson Noise Filter fh{xo). The input values of the 
algorithm are Y (x) , x € 1 (the image) and two numbers m = (27V + 1) x 
{2Nri + 1) and M = {2Nh + 1) x {2Nh +1). In order to improve the resuhs 
we introduce a smoothed version of the estimated similarity distance 



where 



Pk,.o(^)=|JE ^iy)\{Y{y)-Y{Ty)\'-^2f{xo)\ , (27) 
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j^,(,,/z)^exp(^ ^"!,r/""M , (29) 



As smoothing kernels K{y) we can use the Gaussian kernel 

y- 

the following kernel: for y € \Jxo,ri, 

k—uiax{l.j) 

if \\y — xolloo — for some j G {0, 1, ■ • • , Nr}}, and the rectangular kernel 

KAy)^\r^'\l^"''' (31) 

|_ 0, otherwise. 

The best numerical results are obtained using K{y) = Koly) in the definition 
of Pk,xo- Also note that throughout the paper, we symmetrize the image near 
the frontier. 

We present below the denoising algorithm which realizes OWPNF and 
shows its performance on some test images. 



Algorithm: Optimal Weights Poisson Noise Filter (OWPNF) 

First step: 

Repeat for each € I 

Let a = 1 (give the initial value of a) 
compute pK,a:o (2^0 by 
reorder pK,xo{xi) as increasing sequence 
loop from fc = 1 to M 

/(^o)+Et.pU(--) fa,) thong - 

else quit loop 
else continue loop 
end loop 

compute w{x.) = 

compute /'(xo) = Y^xi^v^^.u w{xi)Y{xi). 
Second step: 

For each xq G I, compute 7(2:0) = Y.xe\J^„ k ^' 
If7(xo)<5 

compute 7(.o) = 

else /(a:o) = /'(xq). 



12 



Qiyu JIN et al. 




NM1SE=0.0739 
M = 19 X 19 
m = 13 X 13 
(a) Spots 



NMISE=0.0618 
M = 15 X 15 
m = 5 X 5 
(b) Galaxy 



NM1SE=0.0368 
A/ = 9 X 9 
m = 19 X 19 
(c) Ridges 



NMISE=0.1061 
M = 15 X 15 
m = 21 X 21 
(d) Barbara 



NMISE=0.0855 
M = 11x11 
m = 17 X 17 
(e) Cells 



Fig. 2 These images restored by the first step of our algorithm. 



Algorithm 


Our 


EUI 


MS- VST 


MS- VST 


PH-HMT 


Algorithm 


algorithm 


+BM3D 


+ 7/9 


-1- B3 




Spots[0.08,4.99] 


0.0259 


0.0358 


0.0602 


0.0810 


0.048 


Galaxy[0, 5] 


0.0285 


0.0297 


0.0357 


0.0338 


0.030 


Ridges[0.05,0.85] 


0.0162 


0.0121 


0.0193 


0.0416 




Barbara[0.93, 15.73] 


0.1061 


0.0863 


0.2391 


0.3777 


0.159 


Cells[0.53, 16.93] 


0.0794 


0.0643 


0.0909 


0.1487 


0.082 



Table 2 Comparison between Optimal Weights Filter, MS- VST + 7/9, and MS- VST -|- 
B3. 



Note that tlie presented algorithm is divided into two steps: in the first step 
we reconstruct the image by OWPNF from noisy data; in the second step, we 
smooth the image by a Gaussian kernel. This is explained by the fact that 
images with brightness between and 255 (like Barbara) are well denoised 
by the first step, but for the low count levels images, the restored images by 
OWPNF are not smooth enough (see Figure [5]). For these types of images, 
we introduce an additional smoothing using a Gaussian kernel (see the second 
step of the algorithm). 

Our numerical experiments are done in the same way as in (ssj and plj 
to produce comparable results; we also use the same set of test images (all 
of 256 X 256 in size): Spots [0.08,4.99], Galaxy [0,5], Ridges [0.05, 0.85], Bar- 
bara [0.93,15.73], and Cells [0.53,16.93]. The authors of (H and ^ kindly 
provided us with their programs and the test images. 

Figures[3]-[7]iUustrate the visual quality of the denoised images using OW- 



PNF,EUH-B]VI3D [23, IMS- VST 7/9 [33|, MS- VST -h B3 [331 and PH-HJVIP 
0. 

Table [2] shows the NJVIISE values of images reconstructed by OWPNF, 
EUI+BJVI3D, MS- VST + 7/9, and MS- VST + B3. For Spots [0.08,4.99] and 
Galaxy [0, 5], our results are the best; for Ridges [0.05, 0.85], Barbara [0.93, 15.73], 
and Cells [0.53, 16.93], the method EUI+BM3D gives the best resuhs, but our 
method is also very competitive. 
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4 Proofs of the main results 

4.1 Proof of Theorem □ 

We begin with some prehminary resuhs. The fohowing lemma is similar to 
Theorem 1 of Sacks and Ylvisaker [s^ where, however, the inequality con- 
straints are absent. 

Lemma 1 Let gp{w) he defined by ill]) . Then there are unique weights Wp 
which minimize gp{w) subject to given by 

1 



Wp{x) 



{b-Xp{x))+, 



where b and A are uniquely determined by the following two equations: 

E -j^{b-Xp{x)y^l, 

E 7^(6-Ap(x))+p(a:) = A. 

Proof Consider the Lagrange function 



(32) 

(33) 
(34) 



G{wM,h)^gp{w)-2h\ E w{x)-l\-2 ^ Kx)w{x), 



Q.h 



where &o G K and b e R'^^^'^i^^a.h) jg ^ vector with components b{x) > 0, 
X G Vxo.h- Let Wp be a minimizcr of gp{w) under the constraints ([7]). By 
standard results (cf. Theorem 2.2 of Rockafellar (1993 [1^); see also Theorem 
3.9 of Whittle (1971 [12)), there are Lagrange multipliers 6o e K and b{x) > 0, 
X e Vxg^h such that the following Karush-Kuhn- Tucker conditions hold: for 
any x G U^o,/i, 



d 



dw (x) 



Giw) 



with 



and 



= 2Xp{x) + 2f{x)wp{x) -2b- 2b{x) = 0, 



d 



db {x) 



G{w) 



E Wp{y)-1 = Q, 

0, if6(x)>0, 



Wp{x) 



>0, \ib{x) =0. 



(35) 
(36) 

(37) 
(38) 
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(Notice that the gradients of the equality constraint function h (w) — X^xeu ^ 
w{x) — 1 and of the active inequahty constraints (w) = w (x) , x € 'Vxo,h, 
are always linearly independent, since the number of inactive inequality con- 
straints is strictly less than card Vxo.h-) 

If b (x) = 0, then by (j38p we have Wp (x) > 0, so that by ([55]) we obtain 
b — \p{x) = f{x)wp{x) > and 

ib-Xp{x)) + 

Wp{x) = 



If b {x) > 0, then by ([57| we have Wp (x) = 0. Taking into account we 
obtain 

b - Xp{x) = -b{x) < 0, (39) 

so that 

(6-Ap(x))+ 



Wp{x) = 



As to conditions ([55]) and (IMl) . they follow immediately from the constraint 
([57)1 and the equation 

Since the system and has a unique solution (this can be verified 
by substituting j = a), the minimizer of gp {w) is also unique. The assertion 
of the Lemma is proved. 

Now we turn to the proof of Theorem[TJ Applying Lemma[T]with a = b/\, 
we see that the unique optimal weights Wp minimizing gp{w) subject to (O, 
are given by 

wp^j^^{a~p{x))+. (40) 

Since the function 



is strictly increasing and continuous with Mp(0) = and lim Mp{a) = +cxi, 
the equation 

Mp{a) = 1 

has a unique solution on (0,00). The equation (j34p together with (PO)) imply 



A New Poisson Noise Filter based on Weights Optimization 



15 



4.2 Proof of Theorem [2] 

First assume that p (x) = pf_xa{x) — 1/ {x) — f {xo) \ ■ Recall that gp and Wp 
were defined by pT|) and HH). Using the Holder condition (I^Ul) we have, for 
any w, 

gp{wp) < gp{w) < g{w), 
where ^ 

g{w)=l ^ w{x)L\\x ^ xoWL] +^ E ^'(^) (41) 

and _r is a constant satisfying (1211) . Denote w = argminu,^(w). Since Wp 
minimize gp{w) and p{x) < L||a; — a;o||^, we get 

gpiwp) < gp{w) < g{w). 

By Theorem [H 

w{x)^{a-L\\x^x,fS / E {a-L\\x' -xofS . (42) 

where a > is the unique solution on (0, oo) of the equation Mh{a) = 1, where 
Mh{t)= E YL\\x-x4t{t-L\\x-xo\\l)+ >Q. (43) 

Now Theorem [2] is a consequence of the following lemma. 

Lemma 2 Assume that p{x) — L\\x — xq\\^ and that h > cqu^"' with < 
a < orh = Con" uiith cq > ci {L, (3) = (^iM+m±llL^ ^ . Then 

a = C3n-P/(^^+^\l + o{l)) (44) 

and 

g(u') < C4n 2+2^(1 + 0(1)), (45) 
where C3 anc? C4 are constants depending only on (3 and L. 

Proof We first prove (jH]) in the case where h — 1 i.e. Vxo.h = I- In this case 
by the definition of a, we have 

Ml (a) = J2 - - ^o|li)+i||x - xolli = 1. (46) 

Let h = {^jLfl^. Then a - L||x - xqH^ > if and only if \\x - xo|| < h. So 
from we get 

L^T^^ ^ ^ ||x-xo||L^=r. (47) 
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By the definition of the neighborhood U^^ j^, it is easily seen that 



Nh -r/3+2 

h 



and 



Nh 1-20+2 

h 



y: II- - -ofi = sN-'^ Y k"'^" - ^^'^ (1 + ° (1)) • 

l|a;-2;olloc<'i 



Therefore, (|T7| imphes 



7- — -N^h (l + o(l)) = r, 

(/3 + 2) (2/? + 2) ^ ^ 

from which wc infer that 

7^ = cin"W2(l + o(l)) (48) 



with ci = ( r '''^'^g^^a^"^^'' J ^''^^ ■ From and the definition of /i, we obtain 



a = Lh= Lcfn^W (1 + o(l)), 

which proves (|3H) in the case where h = 1. 

We next prove which imphes (HH), under the conditions of the lemma. 
First, notice that iih<h<l, then M;i(i) = Mi_(t), > 0. If ft. > con"", 
where < a < 213+2 ' then it is clear that h > h, for n sufficiently large. 
Therefore Mh (a) = Mi (a), thus we arrive at the equation from which 
we deduce If ft > cqm" and cq > ci, then again ft > ft for n sufficiently 
large. Therefore, (o) = Mi (o), and we arrive again at pSj) . 

We finally prove P5)) . Denote for brevity 

^ (ft^-||:r-:ro||^)+. 

|a: — 2;o II 00 

Since ft > ft, for rt sufficiently large, we have (a) = Af (a) = 1 and 
Gh = Therefore by gl]), gl]) and (02]), we get 



M^(a) + E||._.„|u<x((«-i|l---o||^)"] ^ 

g{w) = r 



L2G| L 
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Since 



lla; — a^o II oo'^h 



k<Nh l<k<Nh 

-N^h'^\l + o{l)) 



/3 + 2 

1^ ^2-(/3+2)//3 Q , 



we obtain 



5 (^) = r^-^i^L'"''^ (1 + ° (1)) ^ ^4"-^ (1 + o (1)) , 

where C4 is a constant depending on /3 and L. 

Proof of Theorem [2j As p (a;) = \ f {x) — f {xo) \ + 5„, we have 

W{x)p{x)\ <j ^i^)\fi^)-fM\+Sn\ 

<2 Y wix)\fix)-fixo)\\ +251 

Since /(x) < -T, with gp and g by (fTT|) and (|4T]) . we get 

5pH <2ffH+252. 

So 

ffpK)<.9p(ii')<2g(W)+252. 
Therefore, by Lemma [5] and the condition that Sn — O (^n^ 2/+2 j ^ -^ve obtain 

This together with ^ give (|^ . 
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4.3 Proof of Theorem 12] 



Let M' = cardU;^_^ = {2Nh + 1)2/2, m' = cardU^'^^^ = {2Nj] + 1)2/2. 
Denote A^„^^ (y) = f{y) - f[Ty) and r/ (y) = e{y) - e{Ty), where e is defined 
by Q- It is easy to see that 

^ ^ iY{y)~YiTy) f = ^ ^ (y) + i;5(:r) + /(.To) + /(x), 



veu" 



where S'(x) = 52(2;) — Si{x), with 

5i(x)=2 Zi,„,,(y)ry(y), 



Denote 



Then 



52(x)= {rfiy)-fixo)-f{x)) 



l^^i; ^ Zi2 + 



< 



V + fix„) + f{x)-J2f{x„) 



Using one-term Taylor expansion, we obtain 



< 



< 



V + fixo) + fix)-'^2f{xo) 

{fM + f{x) + ev)y^ 

(7(a;o)+7(a;)+ ^1^)1/2 



fixo) + fix)^^^2fixo) 

7{x)-l{xa) 



f{xo) + f{x) + j2fixo) 



Since f{x) > 1/lnn, x G I, (03) and 1^ imply that 
^ I^h^^+^^Six)[ 



(49) 



(50) 



(51) 



We shall use three lemmas to finish the Proof of Theorem [31 
The following lemma can be deduced form the results in Borovkov Q, see 
also Merlevede, Peligrad and Rio [23| . 
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Lemma 3 //, for some 6 > 0,j £ (0, 1) and K > 1 we have 

supEexp((5|X,|'') <K,i = l,...,n, 

then there are two positive constants ci and C2 depending only on S, 7 and K 
such that, for any t > 0, 

P [ ^ X, > t j < exp (-cii^/n) + nexp (-cat^) . 



Lemma 4 Assume that h = cqu ^1^+^ with cq > ci = (^F ''^"'"§^^2^"'"^^ ^ ^'^ 
_ 1 

and that rj ~ C2n 2/3+2 _ Suppose that the function f satisfies the local Holder 
condition I120\} . Then there exists a constant C4 > depending only on /3 and 
L, such that 



'I max Pxoix) > CiU '^Vlan^^o(^n ^1+2^ 



(52) 



Proof Note that 



Ee^fe) ^ ^e'=:L±i^ = ef{y)e^'~'^f^y^ < ere^^-^)^. 



kl 

k=0 



From this inequahty we easily deduce that 

supE(el^(^)l''') < sup (Ee^(^--^)+^(^)+V2(7^ 

y \ ^ y ^ 



where 

Z{y) = 2Zi,„,, (y) 77 (y) + (ry^ {y) _ /(a^o) - Hx)) . 

By Lemma 131 we infer that there exists two positive constants C5 and cg such 
that 

P { — \S{x)\ > z/V^] < exp(-C5Z^) + m' exp(-cJV^z)i). (53) 
\m' J 



Substituting z ~ ^^^XvLrn'M' into the inequahty (j53p . we see that for m! 
large enough, 



— \S{x)\ > Xfi^ ] < 2 exp (- In m'M') 

m' Jm' ] m'M' 



2 
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From this inequality we easily deduce that 



max \S(x)\ > 



a;eU' , m' 





In m'M' 






\/m' 




Inm'Af 







to' 



Taking M' = (27V/i + 1)2/2 = cgn W /2 and to' = {2Nt]+ 1)^/2 = c^nW2/2, 
we arrive at 

P(B)<C7n~W2, (54) 

where B = {max^jgu' ^ ;;7 |'5'(a;)| > cgn 23+2 ^Inn} and cs is a constant 
depending only on /? and L. Since on the set B we have 

2 „,\^/' 1 



, < -= (55) 

in n J Jin n 



for n large enough, combining (ISTj) . ([M]) and ([55|) . we get ([5^. 
Lemma 5 Suppose that the conditions of Theorem\^ are satisfied. Then 
F (nifhixo) - fixo)\''\Y{x),xe I^'J > C9n-W2 inn) = ©(n-^lrs), 
where cg > is a constant depending only on (3 and L. 

Proof Taking into account ([^ . ([M)) and the independence of e(x), we have 
IE{|A(xo)-/(a;o)nr(x),xelL'„} 



2 



< ^ w{x)p{x) +f{x) ^'(^ 

Since p(a;) < Lh^ , from ([55)) we get 

E{|A(xo)-/(xo)nr(x),xGlL'J 

< j ^ ^;(a;)L/i^ j + /(x) J] «}2(x) 

2 



(56) 



/ 



< 



A New Poisson Noise Filter based on Weights Optimization 



21 



Recall that iju{x) stand for the optimal weights, defined by (psj) . Therefore 



E{\Mxo)-f{xo)\'\Y{x),xeI'^ 
( ( 

E w;i(a;)p:roO 



< 

where vj\ — argniin^]^(w) with 



(57) 



r E ^'(^)- 



Since by Lemma 21 



max pj:g{x) < c^n ■2?+^ V\nn > = 1 ~ 0{n 2/3+2)^ 

the inequality (1571) becomes 
P (E{\fh{xo) - /(xo)Hr(x),x e I'iJ < g,{wi) + 2cln-^^ Inn + L^h^^ 

= 1 - 0(71 2??+2). 

Now, the assertion of the theorem is obtained easily if we note that /i^^ — 

2Q 2/3 _ 2^ 

Cj^qTi 2/3+2 ^-^(Ti;^) < Clin 2/3+2, for some constant C12 depending only on 
/? and L (by Lemma 121 with U^^ ^ instead of \]xo,h)- 

Proof of Theorem ini Using dSH), the condition ^ and bound /(a;) < T 
we obtain 



|A(a:o) - f{x^)\''\Y{x),xe I" , < gi(u;) < ci2. 



for a constant C14 > depending only on /3, L and -T. Applying Lemma O we 
have 



E 



{\fh{x^)~f{x^)?\Y{. 



x),x G I'' 



<P (£'{|/h(a;o) - /(a;o)ni^(a;),a; e 11} < ""+2 In n ) cg?^ 2^+2 inn 



2Ji_ 



+ ^[E{\Uxo) - f{x^)\'\Y{x),xe I" } > CgTi-^ Inn ci 



=0 { n 2fi+2 In 



n , 



where the constant in O depending only on /3, L and Taking expectation 
proves Theorem 121 
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(a) ()ri.t;;inal image (h) iXoisy image (e) OW i-' 




(d) EUI+BM3D (e) MS- VST + 7/9 (f) MS- VST -|- B3 



Fig. 3 Dcnoising an image of simulated spots of different radii (image size; 256 X 256). (a) 
simulated sources (amplitudes G [0.08,4.99]; background = 0.03); (b) observed counts; (c) 
Optimal Weights Filter {M = 19 X 19, m = 13 X 13, d = 2 and H = 1, NMISE = 0.0259); 
(d) Exact unbiased inverse + BM3D {NMISE = 0.0358) (e) MS- VST -I- 7/9 biorthogonal 
wavelet (J = 5, FPR = 0.01, N^ax = 5 iterations, NMISE = 0.0602); (f) MS- VST + B3 
isotropic wavelet (J = 5, FPR = 0.01, Nmax = 5 iterations, NMISE = 0.81). 




(d) EUH-BM3D (e) MS- VST + 7/9 (f) MS- VST + B3 



Fig. 4 Denoising a galaxy image (image size: 256 X 256) . (a) galaxy image (intensity G [0,5]); 
(b) observed counts; (c) Optimal Weights Filter (M = 15 X 15, m = 5x5, (i = 2 and 
H = I, NMISE = 0.0285); (d) Exact unbiased inverse + BM3D {NMISE = 0.0297) 
(e) MS- VST -I- 7/9 biorthogonal wavelet (J = 5, FPR = 0.0001, Nmax = 5 iterations, 
NMISE = 0.0357); (f) MS-VST + B3 isotropic wavelet (J = 3, FPR = 0.0001, Nmax = 10 
iterations, NMISE = 0.0338). 
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(d) EUI+BM3D (e) MS- VST + 7/9 (f) MS- VST + B3 

Fig. 5 Poisson dcnoising of smooth ridges (image size: 256 X 256). (a) intensity image (the 
peak intensities of the 9 vertical ridges vary progressively from 0.1 to 0.5; the inclined ridge 
lias a maximum intensity of 0.3; background = 0.05); (b) Poisson noisy image; (c) Optimal 
Weights Filter (M = 9 X 9, m = 19 X 19, d = 3 and = 2, NMISE = 0.0162); (d) Exact 
unbiased inverse + BM3D {NMISE = 0.0121); (e) MS- VST + 7/9 biorthogonal wavelet 
(J = 5, FPR = 0.001, Nmax = 5 iterations, NMISE = 0.0193); (f) MS- VST -|- B3 isotropic 
wavelet (J = 3, FPR = 0.00001, Nmax = 10 iterations, NMISE = 0.0416). 




(d) EUI-I-BM3D (e) MS- VST -|- 7/9 (f) MS- VST -I- B3 



Fig. 6 Poisson denoising of the Barbara image (image size: 256 X 256). (a) intensity image 
(intensity £ [0.93, 15.73]); (b) Poisson noisy image; (c) Optimal Weights Filter (M = 15 X 
15, m = 21 X 21 and d = 0, NMISE = 0.1061); (d) Exact unbiased inverse -I- BM3D 
{NMISE = 0.0863) (e) MS- VST -|- 7/9 biorthogonal wavelet (J = 4, FPR = 0.001, 
Nmax = 5 iterations, NMISE = 0.2391); (f) MS- VST -|- B3 isotropic wavelet (J = 5, 
FPR = 0.001, Nmax = 5 iterations, NMISE = 0.3777). 
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(d) EUI+BM3D (e) MS- VST + 7/9 (f) MS- VST -|- B3 



Fig. 7 Poissoii dcnoising of fluorescent tubules (image size: 256 X 256). (a) intensity image 
(intensity £ [0.53, 16.93]); (b) Poisson noisy image; (c) Optimal Weights Filter (M = 11 X 11, 
m = 17 X 17, d = 1 and H = 0.6, NMISE = 0.0794); (d) Exact unbiased inverse + 
BM3D {NMISE = 0.0643) (e) MS- VST + 7/9 biorthogonal wavelet (J = 5, FPR = 
0.0001,Nmax = 5 iterations, NMISE = 0.0909); (f) MS- VST + B3 isotropic wavelet (J = 5, 
FPR = 0.001, Nmax = 10 iterations, NMISE = 0.1487). 



